import math
from fractions import Fraction
import numpy as np

np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})#设置np.array用分数表示
p01= np.array([1/3,1/3,1/3]) #状态空间概率
# p= np.array([[1/2,1/2], [1/3,2/3]])
p= np.array([[3/4,1/4,0], [1/4,1/2,1/4],[0,3/4,1/4]])
# p= np.array([[0.5,0.5,0,0], [1,0,0,0],[0,1/3,2/3,0],[1/2,0,1/2,0]])#转移概率矩阵(Transition Probability Matrix)#天气预报下雨问题
# p= np.array([[1/2,1/4,1/4], [2/3,0,1/3],[3/5,2/5,0]])#转移概率矩阵(Transition Probability Matrix)
n=4#转移矩阵步数

#计算p^n
#计算n步转移概率矩阵-p^n(考虑重复的自环)(条件概率)
#n步转移概率矩阵对应条件概率-最大下标对应0+n，条件确定行
def trans_matrix(p,n):
    p_tmp=p
    print("p({})=\n{}".format(1, p_tmp))
    for i in range(n-1):
        p=np.dot(p,p_tmp)
        print("p({})=\n{}".format(i+2,p))
    return p
trans_matrix(p,n)

#s----计算n步首达概率-fij^n(考虑重复的自环)，n=f的行数
#最大下标对应0+n，条件确定行
#返回值：1-n步首达概率
def trans_matrix_noloop(p):
    s=[]
    n=p.shape[0]
    p_tmp=p
    s.append(p_tmp)
    for i in range(n-1):
        p = np.dot(p, p_tmp)
        s.append(p)

    for i in range(n - 1):
        flag=0
        deta = 0#待减去差值，第二个元素要减去第一个的平方（自环的影响）
        tag=1##后续概率未0标志；
        for k in range(len(s)):
            if k==0:
                # s[k][i][i]=s[k][i][i]
                deta = s[0][i][i] * s[0][i][i]

            elif s[k][i][i]!=flag and tag:
                # print("s[{}][{}][{}]>=flag:即{}>={}".format(k, i, i, s[k][i][i], flag))
                flag = s[k][i][i]
                # print("flag={}".format(flag))
                # print("deta=",deta)
                s[k][i][i]=abs(s[k][i][i]-deta)
                deta=s[k][i][i]
                # print("s[{}][{}][{}]={}".format(k,i,i,s[k][i][i]))
                if s[k][i][i]==s[k-1][i][i] or s[k][i][i]==0:##后续概率为0的情况；
                    # print("s[{}][{}][{}]==s[{}][{}][{}],{}={}".format(k,i,i,k-1,i,i,s[k][i][i],s[k-1][i][i]))
                    # s[k+1][i][i]=0
                    tag=0
                # if s[k][i][i]==0 or s[k][i][i]==s[k-1][i][i]:
            #
            else :
                #     # s[k][i][i]==s[k-1][i][i]:#后续概率未0的情况；
                # print("[k][i][i]==s[k-1][i][i]",[k][i][i],s[k-1][i][i])
                # # flag = s[k-1][i][i]
                s[k][i][i]=0
                # flag=float("inf")
                # print("s[{}][{}][{}]={}".format(k,i,i,s[k][i][i]))
    return s

x=np.array(trans_matrix_noloop(p))#1-n步的所有转移概率矩阵
print("1-",x.shape[0],"步首达概率矩阵分别为：\n",x,"形状为：",x.shape)
x=x.sum(axis=0)
print(x.shape[0],"步首达概率之和（即周期fii）为（看对角元素即可）：\n"
        "    (1.遍历态需满足周期为1)\n"
        "    (2.如果1-n步的值均为等比数列，那么其n步的和也为1)\n"
        "    (3.正常返判断即满足u小于无穷的：一般看1-n步对角上的值逐渐变小即满足)\n",x)#
#e----计算n步首达概率-fij^n(考虑重复的自环)

#s---求马尔可夫链的平稳分布及各状态的平均返回时间
import numpy as np
from scipy.linalg import eig
transition_mat = np.matrix([
    [0.7,0.1,0.2],\
    [0.1,0.8,0.1],\
    [0.05,0.05,0.9],\
    ])

S, U = eig(transition_mat.T)
stationary = np.array(U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)
stationary = stationary / np.sum(stationary)
print("平稳分布为：",stationary)
print("各状态的平均返回时间为：",1/stationary)
#e---求马尔可夫链的平稳分布及各状态的平均返回时间